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Introduction 


The  overall  motivation  and  goal  behind  the  current  project  was  to  understand  the  thennal 
transport  in  disordered  networks  in  disordered  networks  at  the  molecular  level  using  molecular 
dynamics  simulations  and  connect  it  to  microscopic  level.  The  next  likely  step  would  be  to  take 
this  understanding  to  Finite  Element  Modeling  level  so  that  this  could  be  applied  at  macroscopic 
epoxy  models  to  study  their  thennal  behavior.  This  approach  will  result  in  making  better 
thennally  efficient  epoxy  composites  with  pre-determined  thermal  and  mechanical 
characteristics. 

Due  to  the  inability  of  pursuing  the  simulations  because  of  unavailability  of  computer  resources, 
during  the  first  three  months  of  the  project  (Oct  2006-Dec  2006),  research  had  been  mostly 
focused  on  reviewing  and  understanding  published  literature  regarding  various  thennal 
conductivity  and  multi-scale  modeling  techniques.  “Simulations  of  epoxy  resins/networks”, 
“Simulations  of  thermal  conductivity  of  ordered  systems”,  “Multi-scale  modeling  of  networks”, 
“Thermal  conductivity  simulations  in  disordered  systems”. 

From  the  literature  search,  we  discovered  that  the  subject  of  understanding  thermal  conductivity 
in  disordered  systems  has  rarely  been  touched  from  the  perspective  of  molecular  dynamics 
simulations.  However,  ordered  systems  such  as  “Carbon  Nanotubes”  have  been  investigated  in 
tenns  of  their  thermal  conductivity  and  good  comparison  with  experimental  results  has  been 
found. 

Predominantly,  there  exist  two  procedures  to  calculate  thermal  conductivity,  namely,  a)  Green- 
Kubo  formalism  and  b)  Fourier  law  methodology.  These  methods  will  be  explained  in  detail  later 
on  in  the  report. 

As  mentioned  above,  the  literature  search  was  also  focused  on  reviewing  various  approaches  to 
scale  the  molecular  models  at  coarse-grain  level  to  study  the  systems  at  a  microscopic  level. 
Though,  there  are  many  different  approaches  taken  by  several  authors,  primarily  two  approaches 
have  been  used  to  coarse-grain  the  models  successfully,  namely  a)  Boltzmann  inversion 
technique  by  Kremer  et  al.  and  b)  Force-Matching  algorithm  by  Gregory  A.  Voth.  Although,  the 
literature  search  was  successfully  performed  for  future  simulations,  the  project  mainly  dealt  with 
performing  simulations  atomistic  level. 

Since,  the  thermal  transport  on  the  ordered  systems  has  been  studied  in  literature  in  quite  detail, 
we  first  studied  the  thennal  properties  of  CNTs,  in  order  to  generate  codes  and  familiarize  with 
the  field  of  thermal  conductivity  calculations  from  simulation  perspective. 

After  that,  we  focus  on  thermal  transport  for  disordered  uncross-linked  epoxy  resin  and  curing 
agent-W  (the  constituents  of  epoxy  networks).  Once,  the  simulation  and  thermal  analysis  of  these 
systems  were  convincingly  done,  we  then  focused  on  building  a  crosslinked  network  and  study 
its  thermal  properties  using  equilibrium  as  well  as  non-equilibrium  molecular  dynamics 
simulation  approach.  Later  on,  simulations  of  CNT  with  epoxy  resin  interface  were  also 
performed  to  study  how  thermal  interface  resistance  across  the  CNT/epoxy  resin  boundary. 
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Here,  it  is  worth  mentioning  the  order  of  the  report. 


First  of  all,  the  field  of  heat  conduction  is  discussed  briefly  from  fundamental  and  literature 
perspective.  Then,  molecular  dynamics  simulations  are  introduced  as  an  alternative  approach  to 
experiments  to  tackle  various  molecular  scale  problems  at  nano-  timescales.  Thereafter,  some 
theoretical  background  is  discussed  to  calculate  thermal  conductivity  using  molecular  dynamics 
simulations. 

Once,  the  basic  theory  and  fundamentals  are  presented,  results  from  thermal  conductivity 
simulations  from  CNTs  are  discussed  briefly.  Then,  we  would  switch  gear  to  disordered  systems 
where  we  the  heat  transport  is  studied  in  un-crosslinked  curing  agent  and  un-crosslinked  resin. 

The  report  then  focuses  on  the  crosslinking  procedure  to  build  the  cross! inked  network  or  epoxy 
based  composites.  The  crosslinked  network  is  then  characterized  with  respect  to  various 
thennodynamic  properties.  The  structure  development  during  crosslinking  stage  has  also  been 
focused.  Thereafter,  thermal  transport  in  crosslinked  systems  has  been  discussed  with  the  help  of 
several  analysis  tools  which  will  be  discussed  during  the  later  part  of  the  report. 

Towards  the  end  of  the  report,  the  results  regarding  thennal  interface  resistance  has  been 
discussed  between  crosslinked  epoxy  resin  and  CNT. 


2 


Heat  Conduction 


The  thermal  transport  in  solid  materials  is  predominantly  governed  by  the  phenomenon  of  heat 
conduction.  Heat  conduction  is  a  phenomenon,  governed  by  Fourier  law  and  states  that  the 
amount  of  heat  transporting  through  a  solid  material  per  unit  area  in  a  unit  time  is  directly 
proportional  to  the  temperature  gradient  across  the  material1.  The  proportionality  constant  is 
broadly  known  as  “thermal  conductivity”.  The  inverse  of  thermal  conductivity  is  known  as 
“thermal  resistively”.  With  the  motivation  of  understanding  the  conductive  or  insulative 
properties,  a  great  amount  of  research  has  been  done  to  appreciate  this  phenomenon  in 
significantly  diverse  class  of  materials  which  include  metals,  ceramics,  alloys,  composites 
materials,  etc1  in  past  few  decades. 

In  ordered  materials,  the  heat  is  transported  using  coherent  vibrations,  known  as  phonons".  The 
phonons  are  in  general  divided  into  two  categories  from  the  perspective  of  heat  conduction, 
namely,  acoustic  and  optical  phonons.  It  is  the  acoustic  phonons  which  give  rise  to  heat  transport 
in  such  systems.  These  vibrations  are  low  frequency  vibrations  of  much  higher  wavelength  that 
that  of  IR-spectra.  The  main  concept  that  differentiate  the  acoustic  phonons  which  optical 
phonons  is  that  fact  that  in  former,  the  neighboring  atoms  moves  in  a  same  direction  while  in 
later,  the  neighboring  atoms  moves  in  opposite  direction,  like  in  stretching  vibrations  of  C-H 
bonds. 

On  the  other  hand,  the  thermal  transport  in  disordered/amorphous  system  is  not  discussed  in 
terms  of  “phonon”  as  the  vibrations  get  scattered  quite  easily  before  traveling  any  significant 
distance  due  to  amorphous  nature.  Among  such  disordered  materials,  epoxy  resin  based  thermo¬ 
sets  and  their  composites  presents  a  class  of  materials  which  possess  excellent  mechanical 
properties  and  hence,  are  often  used  in  several  structural  components  in  aerospace  and 
automobile  industry3.  However,  it  is  well  known  from  experimental  results  that  these  materials 
possess  quite  low  thennal  conductivity  and  are  poor  materials  for  thermal  transport4.  The  issue  of 
thennal  transport  through  these  epoxy  networks  becomes  quite  significant  in  situations  where 
they  are  used  as  adhesive  glues  to  join  two  thermally  conductive  materials,  for  example,  in  heat 
exchangers  in  aerospace  applications,  and  hence  presents  a  severe  bottleneck  for  the  thermal 
transport.  It  has  also  been  shown  experimentally  that  dispersing  highly  thermally  conductive 
phases,  such  as  CNTs  and  nano-particles,  does  not  significantly  enhance  their  thermal 
conduction  properties,  even  above  their  percolation  threshold  limit5.  The  main  reason  behind  this 
is  known  as  Kaptiza  resistance6.  It  is  interfacial  thermal  resistance  which  happens  wherever  two 
materials  subjected  to  thennal  transport  are  allowed  to  be  in  contact.  In  such  a  case,  the  phonons 
from  one  material  gets  scattered  at  the  interface,  leading  to  un-efficient  thermal  transport  across 
the  other  material.  Such  situations  primarily  lead  to  two  possible  solutions,  a)  to  develop  a 
thennally  conductive  pathway  within  the  epoxy  network  or  b)  to  tailor  the  intrinsic  thennal 
properties  of  the  adhesive  network.  Here,  the  latter  alternative  presents  a  clear  motivation  of 
understanding  structure-property  relationship  of  epoxy  based  networks,  especially,  in  terms  of 
possible  thennal  transport  and  will  be  studied  as  a  part  of  this  project  from  simulation 
perspective. 
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In  this  regard,  molecular  modeling  provides  a  complimentary  route  to  experiments  to  appreciate 
structure -property  relationship.  Not  only  modeling  helps  in  analyzing  the  system  at  various 
length  scales  in  a  very  cost  effective  manner  without  performing  actual  experiments,  it  also 
provides  foundation  for  developing  accurate  theories  and  in  designing  better  materials  with 
tailored  properties.  Among  various  modeling  techniques,  molecular  dynamics  (MD)  simulations 
(discussed  later  on)7  present  an  opportunity  to  explore  the  system  properties  by  assigning 
interaction  parameters  to  its  constituent  elements  (atoms)  whose  dynamics  is  governed  by 
Newtonian  equations  of  motion.  To  study  epoxy  networks,  several  groups  have  employed 
atomistic8  as  well  as  coarse-grained9  molecular  dynamics  simulations  to  study  their  material 
properties.  However,  the  literature  has  almost  entirely  been  focused  on  studying  the  mechanical 
and  structural  properties  of  the  systems,  presumably  due  to  their  widely  accepted  role  in  building 
mechanically  strong  components  with  no  reference  to  their  thennal  behavior. 

In  fact,  quite  a  number  of  studies  have  been  performed  by  several  research  groups  in  order  to 
appreciate  thermal  properties  at  molecular  level  using  MD  simulations.  Most  of  these  studies 
have  been  perfonned  on  ordered  systems  such  as  Carbon  Nanotubes10,  ionic  salts11,  silicon12, 
metals  ~  etc.  In  such  studies,  thermal  conductivity  is  often  characterized  in  terms  of  propagation 
of  well  separated  acoustic  phonons  through  the  system.  However,  the  thennal  behavior  of 
amorphous  or  disordered  systems  such  as  epoxy  composites  has  scarcely  been  touched  from  the 
perspective  of  atomistic  simulations. 

This  provided  us  with  the  motivation  to  build  a  crosslinked  network  (of  Epoxy  Resins)  and  to 
undergo  one  of  the  first  studies  of  thermal  transport  in  disordered  systems  in  detail. 
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Molecular  Dynamics  Simulations 


Molecular  dynamics  simulation7  is  one  of  several  simulation  techniques  which  is  quite  powerful 
and  gaining  a  lot  of  interest  in  exploring  various  phenomena  occurring  at  pico-  to  nanosecond 
time  scales.  As  these  simulations  consider  atomic  details,  the  effect  of  structure,  molecular 
chemistry,  etc.  can  be  successfully  studied  on  various  physical  properties.  From  the  perspective 
of  cross-linked  networks,  molecular  dynamics  simulations  also  possess  inherent  capability  of 
addressing  interface  issues  in  composite  materials.  These  interfaces  are  relatively  sharp  in  nature 
and  could  only  be  treated  as  an  approximate  effective  intennediate  phase  in  continuum  modeling. 

The  molecular  dynamics  simulations  work  in  the  principle  of  Newtonian  equations  of  motion.  At 
time  scales  of  femto-seconds  (the  basic  timestep  of  MD  simulations),  the  forces  are  calculated 
for  each  atom  based  on  pre-defined  energetic  potential.  This  potential  include  bonded  (bond, 
angle,  dihedral  angle,  improper  angle)  as  well  as  non-bonded  (van  der  Waals  and  electrostatic) 
interactions.  Once,  the  forces  are  calculated  and  current  velocities,  the  new  co-ordinates  of  each 
atom  is  calculated  based  on  basic  Newton’s  equation  of  motion.  There  are  several  algorithms  in 
literature  to  update  the  co-ordinates  and  velocities  of  the  atoms  as  it  evolves  with  time  such  as 
verlet  algorithm,  predictor-corrector  method,  etc. 

In  order  to  appreciate  thermal  transport  from  molecular  simulation  perspective,  MD  simulations 
provide  two  pathways,  also  known  as  equilibrium  MD  (EMD)  and  non-equilibrium  MD 
(NEMD)  simulations.  These  methods  are  based  on  Green-Kubo  formalism14  and  Fourier  law 
formalism15,  respectively.  While  EMD  simulations  are  best  suited  for  studying  homogenous 
systems,  NEMD  simulations  have  an  advantage  when  studying  heterogeneous  systems16, 
although  they  can  also  be  used  to  study  homogenous  systems.  The  theoretical  background  of 
these  systems  is  discussed  next. 
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Theory:  Equilibrium  Molecular  Dynamics  Simulations 


Equilibrium  molecular  dynamics  simulations  are  often  used  to  calculate  several  transport 
properties  using  Green-Kubo  (GK)  approach14.  This  approach  is  based  on  fluctuation-dissipation 
theorem  and  relates  the  fluctuation  properties  of  the  thermodynamic  system  to  its  linear  response 
properties.  In  other  words,  it  provides  a  pathway  to  relate  out  of  equilibrium  properties  (transport 
coefficients)  with  fluctuations  in  equilibrium  properties.  The  transport  coefficients  are  evaluated 
by  integrating  time  autocorrelation  functions  of  microscopic  fluxes  of  equilibrium  properties. 
Few  of  such  examples  include  diffusion  coefficient,  thermal  conductivity,  electrical  conductivity, 
viscosity,  etc.  In  particular,  the  thermal  conductivity,  X,  using  GK  formalism  is  calculated  by 
integrating  time  auto-correlation  function  of  heat  flux  and  is  given  by  following  equation. 

1  00 
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Here  J(t)  is  the  heat  flux  vector  at  time  t.  In  addition,  V  and  T  represent  the  volume  and 
temperature  of  the  system,  respectively,  while  kg  is  the  Boltzmann  constant.  In  tenns  of 
molecular  dynamics  entities,  J(t)  is  written  as 
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Here,  «,•  and  v,  represent  mass  and  velocity  of  atom  i.  u(ry )  is  total  potential  energy  of  atom  i 
whereas  ry  is  the  distance  between  atom  i  and  j.  In  most  cases,  an  analytic  fonn  of  J(/)  is  used  to 
evaluate  J(t)  which  generally  depend  upon  the  form  of  interaction  potential  u(ry)  employed  in  the 
simulations.  In  our  case,  we  have  used  12-6  Lennard  Jones  potential  for  non-bonded  van  der 
Waals  interactions  along  with  Ewald  summation  for  electrostatic  interactions.  In  such  case,  the 
final  expression  for  microscopic  heat  current  vector  J(t)  sum  becomes 
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In  above  equation  F y  represents  the  short  range  van  der  Waals  force  and  real  part  of  Ewald- 
Coulomb  force.  These  forces  are  computed  in  real  space.  On  the  other  hand,  tensor  S  is  evaluated 
in  Fourier  space.  The  elements  of  tensor  S  are  written  as 
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where  a  and  /?  denote  the  direction  in  reciprocal  space,  k  represents  reciprocal  vector,  Rc 
represents  Ewald  parameter  while  i  and  j  are  atom  indices.  In  addition,  Z,  denotes  the  charge  on 
atom  i.  The  detailed  discussion  for  the  inclusion  of  tensor  S  in  heat  current  vector  has  already 
been  discussed  in  literature  previously17.  Here  it  is  sufficient  to  say  that  it  is  introduced  to  avoid 
divergences  arising  from  the  long  range  columbic  interactions  and  takes  care  of  forces  due  to 
long  range  Ewald  interactions  in  Fourier  space18.  The  justification  to  use  relatively  slower  Ewald 
technique  could  be  deduced  from  the  fact  that  such  a  procedure  to  capture  the  long  range 
interactions  efficiently  is  well  documented  in  literature17,18  and  the  references  there  in. 
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Theory:  Non-Equilibrium  Molecular  Dynamics  Simulations 


The  Fourier  law  approach  is  based  on  the  principle  of  heat  conduction  which  states  that  under 
steady  state  conditions,  amount  of  heat  flow  per  unit  area  in  unit  time  is  directly  proportional  to 
the  temperature  gradient  at  the  cross-section.  This  proportionality  constant  is  better  known  as 
thennal  conductivity  and  is  shown  below.  This  method  for  the  calculation  of  thennal 
conductivity  from  molecular  dynamics  simulations  is  also  known  as  direct  method  as  it  is  quite 
analogous  to  experimental  conditions. 


/  dz 


Here,  Q  is  heat  flow  through  the  cross-section,  A  is  the  cross-sectional  area,  At  is  the  time  for 
which  heat  is  flowing  and  dT/dz  is  the  steady  state  temperature  gradient.  Once,  the  heat  flow  and 
temperature  gradient  is  known  from  the  simulation,  the  thermal  conductivity  is  easily  calculated 
using  above  equation.  The  calculation  of  Q  and  dT/dz  from  molecular  dynamics  simulations  is 
discussed  briefly  below. 

First  of  all,  the  system  of  interest  is  build  as  a  long  thin  slab  along  one  direction  and  is 
equilibrated  at  desired  temperature  and  pressure.  Next,  the  central  part  of  the  slab  is  heated  to 
desired  high  temperature,  Thigh  and  is  kept  at  that  temperature  while  the  end  boundaries  are 
cooled  to  desired  low  temperature  Tiow.  In  order  to  keep  the  regions  at  their  specified 
temperatures,  energy  is  continuously  added  and  taken  off  from  hot  and  cold  region,  respectively 
during  the  course  of  the  simulation.  In  doing  so,  a  temperature  gradient  is  established  across  the 
slab.  In  order  to  calculate  the  temperature  gradient,  the  slab  is  divided  into  pre-dcfincd  number  of 
small  slabs  with  equal  thickness.  Thereafter,  the  temperature  of  each  slab  is  calculated  as 
follows19. 


T:  = 


-1 


^NikB  k=\ 


mkvl 


where,  N,  is  number  of  atoms  in  7th  slab.  Furthermore,  calculated  temperature  for  each  slab  7}  is 
averaged  over  the  several  pico-seconds  to  get  a  smooth  temperature  profile.  To  get  the  better 
statistics,  the  temperature  profile  is  further  blocked  averaged  over  several  blocks.  At  last,  the 
temperature  gradient  is  calculated  by  the  slope  of  resulting  temperature  profile. 


On  the  other  hand,  heat  flux  per  unit  area,  Q/AAt  is  calculated  as  follows 


Q  1  /l  V2'  /  2  2s.\ 

-2^mk(yk-vpk) 


AAt  AAt  \  2 


where  vpi,  and  \>k  are  the  velocities  of  the  atoms  before  and  after  rescaling  to  desired  temperature, 
respectively.  Nb  is  the  number  of  atoms  in  the  boundary  layers.  Once  the  temperature  gradient 
and  the  heat  flux  are  known,  the  thermal  conductivity  is  calculated  using  above. 
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Carbon  Nanotube  Simulations:  Approach  and  results  from 
equilibrium  and  non-equilibrium  molecular  dynamics 
simulations. 

First  of  all,  thermal  conductivity  simulations  were  performed  for  (10,10)  armchair  carbon 
Nanotubes  using  EMD  as  well  as  NEMD  simulations.  The  main  reason  to  start  off  with  carbon 
Nanotubes  simulations  initially  was  the  fact  that  its  thermal  behavior  has  been  studied  to  a 
significant  extent  in  literature.  The  nanotubes  were  build  using  Amorphous  cell  builder  in 
Material  Studio.  First  of  all,  the  CNTs  were  equilibrated  using  NVT  and  NPT  simulations. 
Thereafter,  equilibrium  and  non-equilibrium  molecular  dynamics  simulations  were  run  to 
calculate  their  thennal  conductivity.  For  NEMD  simulations,  a  temperature  gradient  was 
established  across  the  Nanotube  while  keeping  track  of  how  much  energy  is  being  put  in  order  to 
maintain  that  temperature  gradient  in  a  steady  state.  The  temperature  gradient  is  shown  in  Figure 
2.  On  the  other  hand,  for  EMD  simulations,  thermal  conductivity  was  calculated  from  the 
integration  of  heat  flux  autocorrelation  function  (shown  in  Figure  3)  as  discussed  in  previous 
section.  The  results  for  both  types  of  simulations  are  listed  as  follows. 

NEMD  approach  (~500W/mK)  @  300K 
EMD  approach  (~1  lOOW/mK)  @  300K 

In  the  literature  the  values  reported  for  the  thermal  conductivity  are  around  -500-3000  W/mK. 
The  spread  in  literature  values  is  attributed  to  method  used,  length  of  the  carbon  Nanotubes 
simulated  as  well  as  the  interaction  potential  incorporated.  The  higher  values  are  mostly  obtained 
with  longer  Nanotubes  with  “Brenner”  interaction  potential  which  considers  short  as  well  as  long 
range  interactions.  It  is  worth  mentioning  that  our  simulations  were  performed  using  “Tersoff’ 
potential  which  only  considers  3-body  short  range  interactions.  The  use  of  Brenner  potential  for 
further  simulations  has  been  left  for  near  future  (due  to  its  unavailability  presently  in  LAMMPS 
molecular  dynamics  simulation  package). 
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Region  Rest  Region  Rest 


Figure  1:  Schematic  of  CNT  for  NEMD  simulations 


(10,10)  armchair  Carbon  Nanotube  (25  nm  length) 

Temperature  Profile 
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(10,10)  armchair  nanotube  (12.5  nm) 

Autocorrelation  Function 
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Figure  3:  Decay  of  heat  flux  autocorrelation  function  for  the  carbon  Nanotube. 


12 


Epoxy  networks  systems  of  interest 


The  system  studied  in  current  study  consisted  of  epoxy  resin,  EPON-862,  curing  agent- W,  also 
known  as  DETDA  and  their  crosslinked  network  and  are  shown  is  shown  in  Figure  4  and  5  along 
with  their  crosslinking  epoxy-amine  reaction.  For  all  simulations  presented  in  this  study,  CVFF 
(Consistent  Valance  Force  Field)  potential  was  used  for  bonded  as  well  as  non-bonded 
interactions.  In  addition,  for  NEMD  simulations,  the  long  range  electrostatic  interactions  were 
modeled  using  PPPM  technique.  However,  since  the  heat  flux  equations  use  Ewald  summation20 
to  calculate  electrostatic  contributions,  the  relatively  slower  Ewald  summation  technique  was 
employed  for  long-range  electrostatic  interactions  for  EMD  simulations.  All  simulations  were 
performed  using  open  source  LAMMPS  (Large-scale  Atomic/Molecular  Massively  Parallel 
Simulator)  software  provided  by  Sandia  National  Labs"  . 

The  initial  configurations  for  both  DETDA  and  EPON-862  were  modeled  using  Material 
Studio®  .  The  system  specifics  after  equilibration  for  each  system  along  with  crosslinked 
network  are  mentioned  in  Table  1.  The  crosslinking  procedure  in  order  to  obtain  a  network  is 
discussed  later  on.  Here,  we  would  like  to  mention  that  we  were  able  to  build  a  crosslinked 
network  with  relaxed  structure  and  was  characterized  successfully  in  terms  of  thermodynamic, 
mechanical  and  structural  properties.  For  current  study,  two  crosslinked  networks  with  ~87.5  % 
crosslinking  was  used  for  both  equilibrium  and  non-equilibrium  simulations. 

Firstly,  all  three  systems  (un-crosslinked  epoxy  resin,  curing  agent  and  crosslinked  network) 
were  minimized  using  conjugate  gradient  method  to  remove  possible  overlaps  between  atoms. 
Thereafter,  a  constant  volume  (NVT)  and  constant  pressure  (NPT)  simulations  were  run  to 
equilibrate  the  temperature  and  density  respectively.  Once  the  initial  equilibrated  system  was 
achieved,  a  series  of  equilibrium  molecular  dynamics  and  non-equilibrium  molecular  dynamics 
simulations  were  perfonned  at  room  temperature. 

For  equilibrium  molecular  dynamics  runs,  after  prior  equilibration,  a  constant  NVE  simulation 
was  perfonned  for  100  pico-seconds  to  get  the  equilibrated  structure  in  micro-canonical 
ensemble.  Once,  the  system  got  equilibrated,  further  constant  energy  (NVE)  simulations  were 
performed  on  newly  equilibrated  system  for  data  collection  for  1.6  ns,  where  data  were  collected 
at  0.01  ps  for  the  calculation  of  heat  current  vector  and  its  subsequent  analysis. 

On  the  non-equilibrium  simulations,  the  elongated  slab  was  divided  into  100  smaller  slabs  of 
equal  thickness.  In  the  central  part  of  the  system,  4  slabs  were  treated  as  hot  region  while  the  2 
slabs  on  each  boundary  were  treated  as  cold  region.  The  system  was  set  to  be  periodic  in  all  3 
dimensions.  For  current  simulations,  the  hot  region  was  kept  at  350K  while  cold  region  was  kept 
at  250  K.  After  the  system  is  equilibrated  from  prior  NPT  simulations,  the  NEMD  simulations  in 
micro-canonical  ensemble  (NVE)  were  run  for  about  Ins  in  order  to  achieve  steady  state  for  heat 
flow.  Once  the  steady-state  was  achieved,  further  simulations  for  1.5  ns  were  run  for  the  data 
collection  at  the  interval  of  0.1  ps  for  further  analysis. 
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EPON-862 


DETDA 

Cyan  (CARBON);  Red  (OXYGEN);  Blue  (Nitrogen);  White  (HYDROGEN) 

Figure  4:  Starting  systems  DETDA  and  EPON-862.  The  depiction  of  colors  is  also  listed 

above. 


Figure  5:  Crosslinking  reaction. 
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Table  1:  Specifications  of  systems  studied 


System  Studied  Equilibrium  MD  Non-Equilibrium  MD 


No.  of 
Atoms 

Dimensions  (A3) 

No.  of  Atoms 

Dimensions  (A3) 

DETDA 

9920 

46.5  x  46.5  x  46.5 

19840 

23.3  x  23.x  373.3 

EPON 

4300 

35.9  x35.9x  35.9 

43000 

35.9x  35.9  x359.9 

7538 

42.56  x  42.56  x 

Crosslinked 

42.56  (90.0) 

15104 

21.4x21.4  X 

System1 

14976 

53.59  x  53.59  x 

337.3  (90.0) 

53.59  (87.5) 

1  For  the  crosslinked  system,  stoichiometric  ratio  of  EPON-862  and  DETDA  were  used.  For  2 
cases  mentioned  in  Table  1,  (128,  64)  and  (256,  128)  molecules  were  used  for  EPON-862  and 
DETDA,  respectively.  The  numbers  in  (  )  signify  crosslinking  %.  After  the  crosslinking 
procedure,  the  un-reacted  entities  were  saturated  with  hydrogens. 
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Results:  Curing  Agent  -  W  (  DETDA)  Thermal  Conductivity 

Equilibrium  Molecular  Dynamics 

Figure  6a)  shows  the  normalized  heat  flux  autocorrelation  function  as  obtained  from  equilibrium 
molecular  dynamics  simulations  in  X,  Y  and  Z  directions.  It  is  observed  from  the  figure  that  the 
autocorrelation  function  decays  down  to  zero  around  ~6-8  ps.  In  addition,  the  figure  also  shows 
that  all  the  autocorrelation  functions  in  different  directions  are  remarkably  similar.  This  could 
lead  to  the  fact  that  the  system  is  indeed  isotropic  and  homogenous. 

Figure  6b)  shows  the  final  thennal  conductivity  in  all  3  different  directions  as  calculated  by 
integrating  the  autocorrelation  function,  where  the  running  average  which  serves  as  a  guide  to 
the  eye  is  shown  in  red.  The  values  of  thermal  conductivity  in  X,  Y  and  Z  directions  were  found 
to  be  0.155  ±0.016,  0.184±0.017  and  0.153±0.009  W/mK  respectively.  The  values  in  all  3 
directions  are  quite  similar  in  nature  suggesting  homogeneous  nature  of  the  system,  providing  the 
average  thennal  conductivity  value  of  0.16  W/mK. 

Non-Equilibrium  Molecular  Dynamics 

Figure  7  shows  the  normalized  velocity  distribution  of  the  carbon  atoms  (Cl,  C2,  and  C3,  as 
discussed  in  previous  section),  in  different  heating  zones.  In  the  figure,  the  velocity  distribution 
in  hot  zone  is  shown  as  black  curve  while  the  grey  curve  shows  the  velocity  distribution  in  the 
middle  of  thennostated  zone.  Both  curves  overlap  each  other  quite  nicely  and  are  in  close 
agreement  with  Maxwell-Boltzmann  law  for  velocity  distribution.  This  similarity,  in  addition  to 
nice  overlap  suggests  that  there  exists  the  local  thermal  equilibrium  within  the  slab  which 
confirms  the  use  of  non-equilibrium  molecular  dynamics  simulations  to  calculate  various  thermal 
properties  as  discussed  by  Chatrenne  et  al.23. 

As  mentioned  before,  to  calculate  the  thermal  conductivity  from  the  NEMD  simulations,  one 
needs  to  make  sure  that  the  system  does  reach  in  a  steady  state  (constant  heat  flux)  before 
estimating  temperature  profile  and  heat  flux.  For  our  simulations,  the  steady  state  was  reached 
after  1.2  ns  as  seen  from  the  similar  slope  of  both  curves  as  shown  in  figure  8. 

The  steady  state  temperature  profile  for  the  DETDA  system  is  shown  in  figure  9  where  the 
temperature  is  plotted  as  the  function  of  layer  number.  Each  layer/slab  was  3. 73 A  was  estimated 
to  be  around  in  thickness.  The  temperature  of  each  such  layer  was  calculated  as  discussed  in 
theoretical  section,  previously.  In  order  to  achieve  better  statistical  results,  the  figure  only  shows 
averaged  half  system  temperature  profile  from  cold  region  to  hot  region.  The  figure  also  shows 
the  calculated  temperature  gradient  as  fitted  using  least  square  fitting  in  the  middle  part  of  the 
temperature  profile  curve.  The  value  of  the  temperature  gradient  is  estimated  as  1.96  K  per  layer 
(or  0.53  K/A).  After  this  evaluation,  the  thermal  conductivity  is  estimated  to  be  0.2  W/mK.  It  is 
quite  remarkable  to  see  that  the  values  from  equilibrium  molecular  dynamics  simulations  and 
non-equilibrium  molecular  dynamics  simulations  are  quite  similar  and  also  serves  as  an  excellent 
guess,  comparing  with  experimental  results. 

Here,  it  is  worth  mentioning  that  although,  the  current  simulations  were  not  focused  on  the  effect 
of  frequency  of  velocity  update  and  the  controlled  (hot  and  cold  region)  slab  thickness  on  the 
thennal  conductivity,  these  variables  should  not  make  any  significant  difference  in  the  values  of 
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predicted  thermal  conductivity.  The  effect  of  such  variables  has  been  discussed  in  literature 
where  the  authors  did  not  find  any  considerable  change  in  the  values  of  thermal  conductivity  as 
the  function  of  slab  thickness  or  frequency  of  velocity  update. 


EPON  Thermal  Conductivity  Results 

The  epoxy  resin  thennal  conductivity  results  will  be  discussed  along  with  crosslinking  network 
later  in  the  report. 
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Figure  6:  a)  Normalized  heat  flux  autocorrelation  function  as  a  function  of  time  for  3 
different  directions.  B)  Calculated  thermal  conductivity  for  3  different  directions.  The  red 


color  represents  the  running  average  over  the  raw  data. 

Figure  7:  Velocity  distribution  from  non-equilibrium  molecular  dynamics  simulations. 
Black:  velocity  distribution  at  the  hot  region.  Brown:  velocity  distribution  at  the  middle  of 
the  thermostated  region. 
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Figure  8:  The  heat  flow  across  the  DETDA  simulation  box  along  the  length. 


Figure  9:  The  averaged  temperature  profile  (over  both  sides)  as  a  function  of  layer 
number.  Black  solid  line:  fit  to  the  points  to  determine  temperature  gradient. 
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Crosslinking  Approach 

Once,  the  thennal  properties  of  the  uncross-linked  systems  were  calculated,  next  step  was  to 
crosslink  the  system.  In  this  regard,  4  different  crosslink  strategies  were  taken  into  account  which 
are  discussed  below.  The  crosslinking  procedure  was  divided  into  3  main  steps,  namely, 
activation,  crosslinking  and  saturation. 

Activation:  In  this  step,  all  potential  reactive  sites  are  activated  for  both  epoxy  resin  and  curing 
agent  as  discussed  by  Xu24.  This  is  done  by  first  hydrating  epoxy  oxygen  atoms  at  the  ends  and 
hence  creating  a  reactive  end  methylene  group.  In  case  of  curing  agent,  DETDA,  both  hydrogen 
atoms  are  removed  to  make  it  chemically  active.  In  addition,  the  partial  charges  are  also  adjusted 
on  these  carbon  and  nitrogen  atoms  in  respective  molecules  due  to  addition  and  removal  of 
hydrogen  atoms  to  keep  the  system  neutral.  The  activation  step  is  shown  schematically  in  Figure 
10. 

Cross-linking  Approaches: 

Approach  1:  This  approach  assumes  equal  reactivity  of  primary  and  second  amine. 

Approach  2:  This  approach  assumes  that  primary  amine  has  much  higher  reactivity  than 
secondary  amine.  In  this  case,  all  the  nitrogen  atoms  are  reacted  once  before  any  secondary 
nitrogen  reaction. 

Approach  3:  In  this  approach,  a  nitrogen  is  selected  randomly  and  is  fully  reacted  (twice)  before 
any  other  nitrogen.  This  approach  was  employed  to  study  dihedral  energy  buildup  during  the 
course  of  the  cross-linking  as  discussed  later. 

Approach  4:  It  is  the  dynamic  crosslinking  approach.  In  this  approach,  all  potential  reactive  pairs 
are  reacted  simultaneously  within  the  cutoff  in  compliance  with  equal  reactivity  assumption. 

The  motivation  for  studying  approaches  2  and  3  in  addition  to  approaches  1  and  4  as  studied  by 
Xu24  and  Heine9,  respectively  was  to  realize  if  they  led  to  better  equilibrated  crosslinked 
structures  with  lower  energies. 

Crosslinking  Algorithm:  The  crosslinking  algorithm  is  shown  in  Figure  11a  and  1  lb  for  cases  1- 
3  and  case  4,  respectively  as  a  flow  chart  and  is  discussed  below. 

Step  1 :  In  this  step,  initial  cutoff  distances  as  well  as  crosslinking  limits  (in  terms  of  %)  are 
defined  for  the  crosslinking  algorithm.  For  cases  1-3,  upper  bound  cutoff  distance  (10  A)  was 
defined  while  for  case  4,  lower  bound  cutoff  distance  (4  A)  was  defined.  Although  the 
crosslinking  limit  was  set  to  be  100%  originally,  the  simulations  could  be  stopped  at  any  time. 

Step  2:  In  this  step,  system  is  equilibrated  using  NPT  simulations  for  40  pico  seconds  (10  ps  of 
equilibration  at  elevated  temperature,  600K;  25  ps  of  annealing  to  300K  followed  by  5ps  of 
equilibration  at  300K).  For  the  detennination  of  equilibration  time,  first  of  all,  RMSD  (root  mean 
square  distance)  of  reactive  carbon  and  nitrogen  atoms  was  estimated.  In  40  pico-seconds, 
RMSD  of  nitrogen  and  carbon  atoms  were  calculated  to  be  ~4.6  and  4.2  A,  respectively.  The 
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estimation  of  equilibration  time  using  this  approach  make  sure  enough  random  movement  of 
reactive  species  before  each  crosslinking  reaction  and  ensures  that  the  reaction  is  indeed 
happening  in  a  pseudo-equilibrium  fashion. 

Step  3:  In  this  step,  all  potential  un-crosslinked  reactive  pairs  are  identified.  In  addition,  all  pairs 
that  pertain  to  possible  catenation  in  the  network  are  removed  from  the  possible  reaction.  For 
cases  1-3,  the  minimum  distance  pair  is  then  selected  and  a  bond  is  formed  if  its  distance  is  less 
than  the  cutoff  distance.  On  the  other  hand,  for  case  4,  all  bonds  within  the  cutoff  distance  are 
created. 

Step  4a:  If  crosslinking  does  occur,  the  topological  infonnation  is  updated  by  introducing  new 
bonds,  angles,  dihedral  angles  and  improper  angles  into  the  system.  Thereafter,  the  topology  is 
relaxed  using  a  multi-step  relaxation  procedure.  This  procedure  essentially  resists  the  abrupt 
change  in  the  co-ordinates  of  atoms  in  the  vicinity  of  newly  fonned  bond  and  ensures  the  slow 
relaxation  of  atoms  comprising  of  newly  formed  topology.  In  this  procedure,  the  force  constant 
of  the  newly  formed  bond  is  gradually  increased  from  zero  to  its  equilibrium  value  in  a  series  of 
steps.  Similarly,  the  equilibrium  bond  length  of  the  new  bond  is  gradually  decreased  from 
original  distance  to  its  actual  equilibrium  value.  The  system  is  minimized  at  each  step  with 
effective  bond  parameters.  For  our  cases,  we  used  five  such  steps  as  shown  in  Table  2. 

Step  4b:  If  crosslink  does  not  occur,  the  system  is  further  equilibrated  according  to  step  2. 
However,  for  case  4,  the  cutoff  distance  is  increased  by  0.25  A  for  the  next  iteration. 

Step  5:  After  the  minimization,  if  crosslinking  limit  is  not  reached,  the  algorithm  jumps  to  step  2 
where  the  equilibration  of  the  new  structure  is  performed;  else  the  simulations  are  stopped. 

Saturation:  Once  the  crosslinking  algorithm  finishes,  the  remaining  active  sites  on  both  un¬ 
reacted  amine  nitrogen  atoms  and  epoxy  carbon  atoms  are  saturated  with  hydrogen.  Here  charges 
on  these  complementary  atoms  are  also  adjusted  to  keep  the  system  neutral.  Although,  it  is  a 
subtle  point,  but  it  is  worth  noticing  that  once  a  crosslink  is  created  between  amine  nitrogen  and 
epoxy  carbon,  the  charge  distribution  around  these  atoms  is  going  to  be  different  than  if  the 
crosslink  was  not  been  formed.  In  order  to  rectify  this  problem,  charges  were  evaluated  on  a 
model  molecule,  as  shown  in  Figure  12  which  has  the  similar  topology  as  around  newly  fonned 
crosslinks.  Then,  the  original  charge  distribution  around  the  new  crosslinks  was  replaced  by 
charge  distribution  shown  in  Figure  12,  in  the  shaded  region  keeping  the  system  neutral. 

In  order  to  achieve  a  relaxed  crosslink  network,  the  system  is  further  equilibrated  at  zero  pressure 
using  constant  pressure  simulations  by  first  heating  at  600  K  for  100  pico  seconds  followed  by 
annealing  to  300K  in  6  nano  seconds. 

Although,  the  above  methodology  has  been  discussed  for  system  sizes  of  (16,  8),  we  also 
performed  crosslinking  simulations  for  relatively  larger  (128,  64)  systems  in  2  different 
geometries  (cube  and  slab)  using  dynamic  crosslinking  approach  as  discussed  above,  the  main 
motivation  being  the  evaluation  of  thermal  properties  as  discussed  later.  The  use  of  this  dynamic 
crosslinking  approach  over  other  approaches  should  become  more  transparent  during  the 
discussion  concerning  comparisons  of  different  approaches  used  to  crosslink  (16,  8)  systems. 
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While,  the  results  from  (16,  8)  simulations  are  used  to  compare  different  approaches  in  terms  of 
buildup  of  intra-  and  inter-molecular  energies  as  a  function  of  crosslinking,  the  results  from  (128, 
64)  simulations  are  used  to  analyze  different  thermal,  thennodynamic  and  structural  properties.  It 
will  also  be  shown  that  some  quantities  such  as  volume  shrinkage  during  crosslinking  are  not 
convincing  extracted  from  (16,  8)  simulations,  however,  they  come  out  quite  nicely  with 
relatively  large  (128,  64)  systems. 

Below,  results  from  various  cross-linking  approaches  are  discussed  for  smaller  (16,  8)  systems. 
In  this  regard,  various  intra-  and  inter-molecular  energies  along  with  their  computational  times 
have  been  listed  in  Table  3  for  un-crosslinked  vs.  crosslinked  system  using  different  approaches. 

The  important  points  from  the  table  are  as  follows 

•  All  energies  except  columbic  energy  are  increased  with  respect  to  un-crosslinked  system. 

•  The  rise  in  bond,  angle  and  dihedral  energies  could  be  attributed  to  newly  formed 
topology. 

•  The  increase  in  van  der  Waals  energy  is  due  to  the  increase  in  number  of  1-4  interactions. 
The  CVFF  force  field  treats  1-4  interactions  as  non-bonded  interactions. 

•  Dynamic  cross-linking  approach  (approach  4)  lead  to  better  equilibrated  crosslinked 
structure  in  terms  of  various  intra-  and  inter-molecular  energies.  Here  almost  all  energies 
are  at  least  among  all  studied  approaches. 

•  Dynamic  crosslinking  approach  also  has  an  advantage  with  respect  to  computational  time 
over  other  approaches.  In  further  support  of  dynamic  crosslinking  approach,  we  would 
like  to  point  out  that  for  the  crosslinking  of  larger  (128,  64)  system  with  256  possible 
crosslinks,  -95%  crosslinking  was  achieved  in  only  75  iterations.  Hence,  it  shows  a 
significant  reduction  in  computational  time  for  crosslinking  (l/3ld)  with  respect  to  other  3 
approaches  which  would  take  more  than  250  iterations  on  the  whole  for  same  amount  of 
crosslinking. 

One  of  the  interesting  things  is  to  observe  the  proximity  of  reactive  sites  during  the  crosslinking 
stage  using  various  approaches.  One  of  such  results  is  shown  in  Figure  13  where  the  initial 
distance  at  which  bond  is  fonned  is  shown  as  the  function  of  crosslinking  %  with  dynamic 
crosslinking  procedure  for  (128,  64)  cube  system.  As  seen  from  the  figure,  about  80%  and  90% 
of  the  crosslinking  has  been  successfully  achieved  within  the  cutoff  distance  of  4  A  and  5.5  A, 
respectively.  When  compared  to  experimental  crosslinking  where  60-70%  crosslinked  is  reached 
generally,  the  figure  clearly  shows  that  using  the  dynamic  crosslinking  approach  as  discussed 
before,  a  well-relaxed  crosslinked  network  can  be  attained  successfully  without  imposing  any 
crosslinking  strain  within  significantly  less  time  as  compared  to  other  approaches.  The  final 
crosslinked  structure  is  shown  in  figure  14. 
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Figure  10.  Activation  of  epoxy  and  amine  ends  for  the  cross-linking  reaction 
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Figure  11.  Crosslinking  algorithm.  A)  for  approach  1-3  and  B)  for  approach  4. 
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Figure  12:  The  shaded  section  of  the  molecule  represents  the  cross-linking  surroundings  in 
which  charges  were  redistributed. 


Figure  13:  Plot  of  initial  bond  formation  distance  as  a  function  of  crosslinking. 
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No  of  steps 

Step  1 
Step  2 
Step  3 
Step  4 
Step  5 


Table  2:  Multi-step  relaxation  procedure 

Effective  Force  Effective  equilibrium 


constant” 


kj  5 


2xkQ/5 


3xkQ/5 


4xk0/5 


bond  length 
r0  +4x(r-r0)/5 

ro  +3x(r  — r0)/5 
ro  +  2x(r  — r0)/5 
r0+lx(r-r0)/5 


Table3:  Comparison  of  various  crosslinking  approaches 


Quantities2 * 4 

Un-crosslinked 

System 

Crosslinked  Systems  (30/32  cross 

dinks) 

Apprch.  1 

Apprch.  2 

Apprch.  3 

Apprch.  4 

Bond  Energy  (per  unit) 

419  (0.44) 

471  (0.48) 

473  (0.48) 

472  (0.48) 

468  (0.48) 

Angle  Energy  (per  unit) 

99  (0.06) 

165  (0.09) 

172  (0.1) 

168  (0.1) 

164  (0.09) 

Dihedral  Energy  (per  unit) 

129  (0.06) 

440  (0.18) 

441  (0.18) 

418  (0.17) 

423  (0.17) 

van  der  Waals  Energy 

580 

774 

758 

774 

754 

Columbic  Energy 

-349 

-355 

-335 

-341 

-350 

Total  Energy 

231 

419 

423 

433 

404 

Computational  Time5 

(iterations) 

N/A 

32 

32 

31 

23 

2 

“&o  defines  equilibrium  force  constant 
r0  defines  equilibrium  bond  length 

4  Energy  units  are  in  Kcal/mol.  per  unit  quantities  for  respective  energies  are  shown  in  (). 

5  The  computational  time  is  given  in  terms  of  number  of  iterations  since  every  iteration  has  same 
equilibration  time  (per  iteration)  irrespective  of  the  approach. 
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Crosslinking  Validations 

For  the  validation,  various  material  properties  were  analyzed  for  (128,  64)  system.  Here,  the 
results  will  be  presented  from  cube  geometry.  However,  it  should  be  mentioned  that  similar 
results  were  found  for  slab  geometry  also. 

One  of  the  properties  which  is  directly  evident  during  the  process  of  network  fonnation  in 
experiments  is  the  volume  shrinkage.  It  has  been  reported  in  literature  that  for  epoxy  based 
networks,  the  change  in  volume  accounts  for  ~7-8%  (shrinkage)  for  densely  crosslinked  systems. 
In  order  to  observe  this  shrinkage,  the  instantaneous  volume  of  the  system  during  the 
crosslinking  simulations  was  tracked  continuously.  In  this  regard,  the  resulted  equilibrated 
volume  as  curing  progresses  is  shown  in  Figure  15.  It  is  quite  interesting  to  note  that  for  the 
crosslinking  of  highly  dense  system  (--90%  crosslinking),  the  volume  shrinkage  of  7%  was 
observed  in  our  simulations  which  is  in  good  agreement  in  experimental  findings,  given  the 
complexity  of  curing  process.  It  is  also  important  to  notice  here  that  such  a  volume  shrinkage 
was  not  convincingly  observed  in  case  of  smaller  (16,  8)  systems  due  to  significant  statistical 
variations  as  shown  in  inset  of  figure  15. 

In  order  to  study  the  temperature  dependence  of  various  thermodynamic  quantities,  an  annealing 
simulation  was  performed  where  the  crosslinked  system  was  slowly  cooled  down  at  the  rate  of 
10K/200ps  from  500K  to  200K  under  atmospheric  pressure  which  was  controlled  by  Nose- 
Hoover  thennostat  and  barostat.  In  this  regard,  the  change  of  density  of  the  cured  network  as  a 
function  of  temperature  is  plotted  in  Figure  16.  The  density  at  room  temperature  is  found  to  be 
1.12  which  is  in  good  agreement  with  experiments  .  Apart  from  the  expected  negative  slope  of 
density  vs.  temperature,  a  kink  is  observed  in  figure  16,  where  the  change  of  slope  occurs.  The 
kink  is  usually  identified  as  the  glass  transition  temperature  (Tg)  which  in  our  case  occurs  around 
105  °C.  The  glass  transition  temperature  is  a  second  order  transition  and  its  value  is  often 
governed  by  cooling  rate  as  well  as  method  of  measurement.  In  case  of  our  systems,  the 
predicted  Tg  from  our  simulations  are  in  fair  agreement  with  experimental  "  as  well  as  simulation 
results  on  similar  systems,  providing  rationalization  for  studied  annealing  rate. 


Coefficient  of  thermal  expansion  is  another  thermodynamic  property  which  is  of  great 
importance  in  composite  materials  as  they  are  used  at  various  temperatures  in  practical 
applications.  The  coefficient  of  volume  thermal  expansion  (CVTE),  a,  is  defined  as 


,dT, 


P 


where,  V,  T  and  P  are  the  volume,  temperature  and  pressure  of  the  system  respectively.  For 
isotropic  materials,  the  coefficient  of  linear  thermal  expansion  (CLTE),  /?,  is  related  to  a  as 


In  this  regard,  the  fractional  change  in  volume  (dV/V)  vs.  temperature  is  plotted  in  Figure  17.  As 
expected,  a  kink  is  observed  around  105  °C,  suggestive  of  glass  transition  temperature. 
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Additionally,  the  slope  from  the  figure  17  can  be  identified  as  CVTE,  a.  From  above  equations, 
the  calculated  values  of  CLTE  were  found  to  be  40  ppm/°C  and  90  below  and  above  the  glass 
transition  temperature.  Given  the  fact,  that  the  annealing  was  done  on  highly  crosslinked  network 
(—90%  crosslinked),  one  could  say  that  the  predicted  values  are  in  fair  agreement  with 
experimental  findings. 

In  addition  to  thennodynamic  properties,  it  is  also  important  to  appreciate  different  kinetic 
processes  during  the  course  of  network  formation.  The  processes  are  governed  by  reaction  rates 
of  various  reactions  and  are  often  analyzed  in  terms  of  molecular  weight  build  up,  its 
distribution,  cycle  ranks,  gel  point,  etc.  using  experimental  techniques  such  as  chromatography, 
etc.  These  processes  eventually  define  the  molecular  structure  of  the  crosslinked  network.  Here, 
we  present  few  of  such  results  in  order  to  support  our  cross-linking  approach  from  reaction 
kinetics  perspective. 

In  this  regard,  weight  average  molecular  weight,  Mw  is  shown  as  a  function  of  cure  conversion  in 
Figure  18.  The  figure  shows  that  most  of  the  rise  in  the  molecular  weight  happens  towards  higher 
cure  conversion  (~0. 65-0. 75).  This  trend  of  molecular  weight  build-up  is  in  accordance  with  the 
principle  of  step-growth  polymerization,  the  same  approach  which  is  used  to  cure  epoxy  resins 
experimentally.  From  theoretical  perspective,  the  figure  also  shows  the  build  of  molecular  weight 
of  largest  group  within  the  system,  the  inception  of  secondary  cycles  and  their  trend,  along  with 
theoretical  gel  point  (0.574)  as  predicted  by  kinetic  theory  for  bi-functional  epoxy  resins  and 
tetra-functional  cross-linkers.  A  secondary  cycle  is  defined  as  an  intra-molecular  reaction 
happening  within  a  group.  The  inception  of  the  secondary  cycles  has  also  been  treated  as  a  gel 
point  in  literature.  From  the  figure,  it  is  quite  clear  that  in  our  simulations,  the  inception  of 
secondary  cycles  occur  around  0.57  providing  nice  agreement  between  simulation  results  and 
theoretical  predictions.  The  gel  point  can  also  be  estimated  from  the  inflection  point  of  largest 
group  molecular  weight  build  up.  Although  it  is  quite  difficult  to  identify  the  inflection  point 
with  great  accuracy  from  figure  18,  one  can  safely  say  that  it  seems  to  be  occurring  between  0.6 
and  0.65,  providing  a  good  estimation  of  gel  point.  Here,  we  would  also  like  to  point  that  that 
towards  higher  cure  conversion  (>0.8),  most  of  the  reactions  happening  are  intra-molecular 
reactions.  This  is  supported  by  the  fact  that  in  this  region,  a  linear  increase  in  number  of 
secondary  cycles  are  observed  as  a  function  of  cure  conversion  along  with  negligible  build  up  of 
molecular  weight. 

It  is  known  from  the  literature  that  near  gel  point,  the  probability  of  reaction  between  largest  and 
second  largest  group  is  maximum.  The  reasons  being,  a)  it  gives  rise  to  sharp  increase  in 
molecular  weight  of  largest  group,  i.e.,  inflection  point  and  b)  after  this  reaction,  most  of  the 
reactions  are  going  to  occur  within  this  largest  group  and  hence  it  also  corroborates  with  the 
concept  of  inception  of  secondary  cycles.  In  this  regard,  Chiu  and  co-workers  performed  Monte 
Carlo  simulations  on  curing  of  epoxides  with  amines  and  estimated  the  gel  point  using  weight 
averaged  reduced  molecular  weight  (RMW)  approach.  RMW  is  defined  as  the  average  molecular 
weight  of  all  reacting  groups  except  the  largest  one.  One  can  easily  foresee  that  the  RMW  first 
increases  with  cure  conversion  and  then  decreases.  Here,  the  maximum  in  RMW  is  also 
characterized  as  a  gel  point.  In  order  to  compare  our  approach  with  Chiu  et  al.,  we  have  also 
estimated  the  gel  point  using  RMW  approach.  The  results  are  presented  in  Figure  19  where 
molecular  weight  of  largest,  second  largest  and  weight  averaged  reduced  molecular  weight  is 
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plotted  as  a  function  of  cure  conversion.  It  is  clear  from  the  figure  that  maximum  in  second 
largest  as  well  as  RMW  occur  within  close  proximity  of  the  theoretical  get  point  which  provided 
us  with  confidence  in  our  cross-linking  approach.  Here,  it  is  worth  mentioning  that  both  peaks  in 
Figure  19  are  quite  broad  due  to  limitations  of  system  size  (only  256  possible  cross-links). 
However,  we  predict  that  the  broadness  of  these  peaks  should  decay  down  with  the  size  of  the 
system. 
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Figure  16:  Plot  of  density  variation  as  a  function  of  temperature 
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Temperature  (  C) 

Figure  17:  Plot  of  fractional  volume  change  as  a  function  of  temperature 


Figure  18:  Plot  weight  average  molecular  weight  (circles),  largest  molecular  weight 
(squares)  and  secondary  cycles  (diamonds)  as  a  function  of  cure  conversion.  The  dashed 
lines  are  guide  to  the  eyes. 
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Figure  19:  Plot  of  largest  molecular  weight  (circles),  second-largest  molecular  weight 
(squares)  and  weight  averaged  reduced  molecular  weight  (diamonds)  as  a  function  of  cure 
conversion.  The  dashed  lines  are  guide  to  the  eye. 
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Reduced  Molecular  Weight 


Thermal  Conductivity  of  Crosslinked  Network  from  Non- 
Equilibrium  Approach 

The  thermal  conductivity  calculations  from  non-equilibrium  MD  simulations  were  done  on  the 
basis  of  Fourier  Law  approach  as  discussed  in  the  previous  section.  It  has  been  suggested  in 
literature  that  for  NEMD  simulations  to  be  justifiable,  there  should  exist  a  local  thermal 
equilibrium  within  the  slab.  For  our  simulations,  this  was  confirmed  by  verifying  the  fact  that  the 
velocity  distribution  indeed  followed  the  Maxwell-Boltzmann  curve  for  thermostated  as  well  as 
un-thermostated  regions  as  discussed  previously.  After  that,  the  heat  flux  and  temperature 
gradient  for  calculated  for  all  studied  systems,  namely,  EPON-86,  DETDA  and  crosslinked 
network.  In  this  regard,  Figure  20  shows  the  input  (energy  put  in  hot  region)  and  output  (energy 
taken  off  from  cold  region)  energies  as  the  function  of  simulation  time  in  a  cumulative  manner 
for  all  3  systems.  From  the  figure,  it  should  be  safe  to  assume  that  the  steady  state  has  been 
reached  for  all  3  systems  around  1 .2  ns.  After  steady  state  system  was  reached,  the  heat  flux  per 
unit  time  across  the  cross-section  was  calculated  from  the  slope  of  each  curve  was  evaluated 
using  least  square  fitting  method  and  is  listed  in  Table  4. 

Along  the  same  lines,  the  resultant  temperature  gradient  due  to  steady  state  heat  flux  is  plotted  in 
Figure  21.  As  mentioned  before,  the  system  was  divided  into  100  slabs  or  layers.  The 
temperature  of  each  layer  was  calculated  using  1000  samples  stored  at  every  0.1  ps  for  100  ps. 
The  temperature  was  further  block  averaged  over  13  such  block-average  temperatures.  From  the 
figure,  it  is  quite  clear  that  the  temperature  change  across  the  slab  in  the  un-thermostated  zone  is 
linear  for  all  3  systems.  The  slope,  also  known  as  temperature  gradient  was  calculated  by  fitting 
the  central  portion  of  un-thermostated  zone  to  avoid  any  boundary  effects  using  least  square 
method.  The  values  of  the  temperature  gradient  for  all  3  studied  systems  are  also  listed  in  Table 
4. 

Thereafter,  the  thermal  conductivity  for  all  three  systems  are  listed  in  Table  2.  The  thennal 
conductivities  of  DETDA,  EPON-862  and  their  crosslinked  network  were  found  to  be  0.20,  0.21 
and  0.30  W/m-K,  respectively.  Although,  the  experimental  values  of  thermal  conductivities  of 
DETDA  and  un-crosslinked  EPON-862  were  not  found  it  literature,  the  thennal  conductivity  of 
crosslinked  epoxy  network  has  been  mentioned  between  0.2  and  0.3  W/m-K.  Considering  the 
fact,  the  calculated  thennal  conductivity  is  indeed  in  good  agreement  with  experimental  findings. 
Hence,  non-equilibrium  molecular  dynamics  simulation  presents  a  positive  route  towards 
studying  and  predicted  thermal  properties  of  disordered  or  amorphous  systems. 
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Figure  20:  Heat  Flux  as  a  function  of  time,  a)  DETDA;  b)  EPON-862;  c)  crosslinked 
network.  The  plots  for  b)  and  c)  are  shifted  by  4  and  8  kcal/mol-A2,  respectively,  for  the 
sake  of  clarity.  For  each  system,  the  above  dataset  represents  the  heat  put  in  at  high 
temperature  region  while  lower  dataset  represent  the  heat  taken  off  from  low  temperature 
region. 


Layer  Number 


Figure  21:  Temperature  profile  across  slab  for  various  studied  systems,  a)  DETDA;  b) 
EPON-862;  c)  crosslinked  network.  The  plots  for  a)  and  c)  are  shifted  by  -15  K  and  15  K 
respectively  for  the  sake  of  clarity. 
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Thermal  Conductivity  of  Crosslinked  Network  from 
Equilibrium  Approach 

Figure  22  shows  the  normalized  heat  flux  autocorrelation  function  (as  discussed  previously)  as  a 
function  of  time  for  the  crosslinked  network  along  with  uncrosslinked  epoxy  resin  and  curing 
agent.  The  figure  also  shows  the  thennal  conductivity  for  all  three  systems  as  evaluated  by 
integrating  the  heat  flux  autocorrelation  function.  It  is  observed  from  the  figure  that  the 
autocorrelation  function  is  oscillating  and  decays  down  to  zero  around  —2-3  picoseconds  for 
EPON-862  and  crosslinked  system  while  it  takes  about  6  pico-seconds  for  DETDA.  The 
difference  in  the  decay  times  could  be  qualitatively  attributed  to  the  rigid  and  planar  nature  of  the 
DETDA  molecule  unlike  EPON  and  the  network  which  contain  flexible  ether  connecting 
linkages. 

The  figure  also  shows  the  trend  of  thermal  conductivity  vs.  time  as  evaluated  from  the  time 
integration  of  heat  flux  autocorrelation  function.  For  all  three  cases,  we  find  that  the  thennal 
conductivity  first  increases,  then  decreases  and  then  becomes  constant  when  the  autoconelation 
function  decays  down  to  zero.  Such  a  trend  in  thennal  conductivity  has  been  previously  observed 
for  oscillatory  decaying  heat  flux  auto-conelation  function28.  From  the  plot,  the  thermal 
conductivities  for  DETDA,  EPON-862  and  their  crosslinked  system  were  found  to  be  0.15,  0.13 
and  0.13  W/m-K  at  studied  conditions. 

Here,  one  can  easily  observe  that  the  predicted  thermal  conductivity  values  using  equilibrium 
MD  are  lower  than  their  non-equilibrium  counterpart.  Here,  one  could  attribute  this  to  the  system 
size  effects.  However,  we  would  like  to  mention  that  we  performed  an  equilibrium  simulation  for 
crosslinked  network  for  system  as  big  as  twice  of  initial  system  with  -15000  atoms  (please  refer 
to  Table  1).  In  that  case  too,  we  find  the  thermal  conductivity  values  to  be  0.13  W/m-K  as  with 
smaller  crosslinked  system.  As  we  employed  the  same  force-field  for  equilibrium  and  non¬ 
equilibrium  simulations,  the  effect  due  to  force  field  doesn’t  seem  to  play  a  significant  role  here. 
Hence,  one  could  argue  that  non-equilibrium  molecular  dynamics  simulations  are  possibly  better 
suited  for  thennal  conductivity  prediction  for  crosslinked  network  systems. 
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Figure  22:  Normalized  heat  flux  autocorrelation  function  and  resultant  thermal 
conductivities  for  various  studied  systems,  a)  crosslinked  network  ;  b)  EPON-862;  c) 
crosslinked  DETDA. 


Table  4:  Results  from  EMD  and  NEMD  simulations  of  studied  systems 


System  Studied 

Heat  Flux 
( kcal/mol-A2-ps ) 

Temperature 

Gradient 

(K/A) 

Thermal  Conductivities  (W/m-K) 

Non-Equilibrium  Equilibrium 
Approach  Approach 

DETDA 

3.08  x  10-3 

0.53 

0.20 

0.15 

EPON-862 

3.08  x  10-3 

0.55 

0.21 

0.13 

Crosslinked 

System 

5.12  x  10-3 

0.60 

0.30 

0.13 
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Power  Spectrum  Analysis 

This  section  presents  the  analysis  of  different  vibrations  in  the  crosslinked  network  in  terms  of 
power  spectrum  of  normalized  velocity  autocorrelation  function  (ACF)  (Fourier  Transform).  The 
power  spectrum  presents  the  partial  density  of  states  for  each  atom  and  provides  a  perceptive 
regarding  the  energy  storage  in  different  parts  of  the  molecule28.  The  partial  density  of  states 

9Q 

(PDOS)  in  terms  of  velocity  autocorrelation  function  can  be  written  as" 

r  Np  I  Np 

DpJj(co)  =  c p(t)cos(cdt)At  where  T pit)  =  2](ui/?(f)«uf/l(0))  /  2](ui/?(0)  •ui/J(0)) 

0  i  /  i 

Here,  Dp J^co)  and  r(oj)  denote  the  PDOS  and  the  normalized  velocity  autocorrelation  function 
for  atom  type  /?.  cp  denotes  the  concentration  of  species  of  type  ft ’  The  total  density  of  states  is 
obtained  by  summing  over  the  partial  density  of  states. 

p 


In  this  regard,  Figure  23  shows  the  power  spectrum  of  velocity  ACF  of  various  atomic  entities  in 
the  crosslinked  network.  For  these  calculations,  simulation  was  run  for  40  picoseconds  where 
data  was  stored  at  each  timestep  of  lfs.  The  power  spectrum  was  evaluated  from  averaged 
velocity  ACF  over  X,  Y  and  Z  directions  as  ACF  was  found  to  be  isotropic  in  nature  for  atomic 
types.  As  one  can  see,  the  figure  presents  several  peaks  for  each  atomic  type.  These  peaks 
represent  different  vibrations  (stretching,  bending,  rocking,  etc.)  in  which  corresponding  atomic 
entity  is  involved.  In  ordered  systems,  these  are  often  termed  as  optical  vibrations  or  phonons 
and  are  indeed  captured  using  spectroscopic  techniques  such  as  IR,  Raman  spectroscopy,  etc. 
These  vibrations  are  known  to  be  rarely  involved  in  thermal  transport  through  the  system.  The 
figure  also  presents  a  low  broad  low  frequency  peak  which  is  present  for  all  atomic  entities.  The 
peak  is  better  depicted  in  figure  24.  Generally,  these  peaks  occur  below  low  frequency  limit  of 
IR  spectroscopy  (-400  cm'1)  and  are  associated  which  acoustic  modes  of  vibrations  which  are 
often  responsible  for  thermal  conduction  across  the  system.  For  our  system,  we  see  one  such 
broad  peak,  present  in  all  power  spectra  of  all  atomic  entities.  The  similar  position  of  the  peak  in 
all  atomic  types  suggests  that  indeed  all  entities  are  coherently  involved  in  such  low  frequency 
vibrations.  However,  the  broadness  of  the  peak  over  several  THz  suggests  a  possible  distribution 
of  relaxation  times  with  different  overlapping  low  frequency  vibrational  modes.  Recently, 
McGaughey  et  al.  did  a  thermal  conductivity  analysis  of  metal  organic  framework  (MOF-5) 
crystal"  .  The  power  spectrum  analysis  from  their  work  give  rise  to  several  sharp  peaks  in  the 
low  frequency  vibrational  regime  (acoustic  region)  unlike  our  work  where  we  observe  only  one 
broad  peak.  Here,  the  differences  in  the  peak  characteristics  in  low  frequency  regime  could  be 
attributed  to  the  disordered  and  amorphous  nature  of  our  system. 
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Figure  23.  Power  spectra  of  velocity  autocorrelation  function  of  various  atomic  entities  of 
crosslinked  network,  a)  sp2  benzene  carbon;  b)  amine  nitrogen;  c)  ether  and  hydroxyl 


't 

oxygen;  and  d)  sp  carbon  in  methyl  and  methylene  groups. 

Figure  24.  Low  frequency  power  spectra  of  velocity  autocorrelation  function  of  various 
atomic  entities  of  crosslinked  network,  a)  sp2  benzene  carbon;  b)  amine  nitrogen;  c)  ether 
and  hydroxyl  oxygen;  and  d)  sp3  carbon  in  methyl  and  methylene  groups. 
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Pair  Contribution  Analysis 


Although  the  results  from  equilibrium  molecular  dynamics  simulations  provide  a  lower  estimate 
of  thermal  conductivity,  it  should  be  interesting  to  study  the  contribution  of  different  interactions 
towards  thennal  conductivity.  The  heat  flux  equation  shows  that  the  heat  flux  indeed  comprises 
of  several  terms.  These  terms  are  associated  with  a)  kinetic  energy;  b)  van  der  Waals  potential 
energy;  c)  electrostatic  potential  energy;  d)  forces  due  to  van  der  Waals  interactions;  e)  forces 
due  to  electrostatic  energy  (within  a  certain  cutoff,  calculated  in  real  space);  and  f)  forces  due  to 
long  range  electrostatic  interactions  (outside  the  cutoff,  calculated  in  Fourier  space). 

In  order  to  capture  different  contributions,  each  tenn  was  evaluated  separately  and  was  stored  for 
each  timestep  during  the  course  of  the  simulation.  As  an  example,  in  order  to  evaluate  thennal 
conductivity  due  to  kinetic  and  potential  energy  contributions  only,  we  summed  up  the  terms 
associated  with  a)  and  b)  when  calculating  heat  flux,  while  ignoring  other  tenns.  For  our 
analysis,  we  represent  the  contribution  or  ignorance  in  terms  of  index  “1”  and  “0”,  i.e.,  “1” 
represents  the  inclusion  of  the  tenn  while  “0”  suggests  that  a  particular  tenn  is  being  ignored 
while  calculating  heat  flux.  For  example,  contribution  description  “100000”  suggests  that  only 
kinetic  energy  contributions  are  taken  into  consideration  while  ignoring  others  during  heat  flux 
contribution.  We  believe  that  such  an  approach  should  take  care  of  possible  cross-terms  as  well 
in  an  implicit  manner. 

In  this  regard,  several  contributions  are  listed  in  Table  5  along  with  the  resultant  thennal 
conductivities  and  are  graphically  shown  in  Figure  25.  On  close  observation,  one  could  conclude 
that 

a)  The  contributions  from  force  tenns  have  a  dominating  effect  than  energy  terms  towards 
final  thermal  conductivity.  Please  compare  rows  (1-4)  vs.  (5-12)  in  Table  5.  Such  a 
behavior  has  also  been  reported  earlier  for  ionic  systems'  . 

b)  The  contribution  from  the  forces  due  to  van  der  Waals  interactions  is  the  most  dominating 
towards  thermal  conductivity  (4th  contribution).  Please  compare  (5-12)  vs.  the  rest. 

c)  The  contributions  due  to  long  range  electrostatic  force  tenns  are  negligible.  Please 
compare  (2,3),  (5,6)  and  (9,10).  This  suggests  that  such  contributions  could  be  neglected 
while  doing  thermal  analysis  for  organic,  neutral  amorphous  systems  with  low  partial 
charges. 

d)  The  contributions  from  the  energy  terms  (kinetic  and  van  der  Waals)  are  low  but 
significant.  Please  see  (1,4).  In  addition,  kinetic  energy  provides  additional  contributions 
while  combining  with  other  terms.  Please  compare  (1+7)  with  10.  This  possibly  suggests 
a  positive  contribution  of  cross-tenns  in  heat  flux  autocorrelation  function  which  involve 
kinetic  energy  contributions. 

Surprisingly,  the  contributions  from  electrostatic  energy  and  short  forces  due  to  short  range 
electrostatic  interactions  were  observed  to  be  negative.  Please  compare  (2  and  4),  (10  and  11) 
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and  (10  and  12).  Such  a  trend  suggests  that  the  cross-terms  involving  the  electrostatic 
interactions  and  forces  are  contributing  negatively  towards  the  heat  flux  autocorrelation  function. 
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Figure  7.  Thermal  conductivity  values  for  different  heat  flux  component  contributions. 
The  details  for  each  contribution  index  are  listed  in  Table  3. 


Table  5:  Thermal  Conductivities  for  different  heat  flux  contributions. 

Contribution  Description  k  (W/m-K)  Standard  Deviation 

1 

100000 

0.0151 

0.0006 

2 

111000 

0.0237 

0.0018 

3 

111001 

0.0238 

0.0018 

4 

110000 

0.0304 

0.0035 

5 

000110 

0.0888 

0.0047 

6 

000111 

0.0889 

0.0047 

7 

011111 

0.1033 

0.0062 

8 

loom 

0.1122 

0.0053 

9 

111110 

0.1313 

0.0071 

10 

linn 

0.1314 

0.0071 

11 

110110 

0.1353 

0.0105 

12 

110100 

0.1392 

0.0083 
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Thermal  Diffusivity  and  Mean  Free  Path  Analysis 


In  addition  to  thermal  conductivity,  thermal  diffusivity  and  heat  capacity  are  two  very  important 
properties  governing  the  thermal  characteristics  of  a  material  and  are  related  through  following 
equation. 


pCP 


where,  D,  X,  p  and  Cp  are  thennal  diffusivity,  thermal  conductivity,  density  and  specific  heat 
capacity  at  constant  pressure.  For  the  estimation  of  thennal  diffusivity,  the  molecular  dynamics 
simulations  provide  a  route  to  calculate  the  heat  capacity  of  a  system  from  the  fluctuations  in  its 
potential  energy  as  shown  below30. 

„  _iu2)~(u)2 


For  our  system,  we  estimated  the  heat  capacity  of  our  system  to  be  1.74  J/g-K  at  300K.  The  heat 
capacity  also  showed  weak  temperature  dependence.  From  previous  estimation  of  thermal 
conductivity  of  0.3  W/m-K  from  non-equilibrium  molecular  dynamics  simulations  along  with 
system  density  of  1.12  gm/cm3,  the  thermal  diffusivity  of  epoxy  networks  was  estimated  to  be 
0.154  mnr/s  at  300K  which  is  in  excellent  agreement  with  experimental  value  of  0.157  mm2/s.27 
The  thennal  diffusivity  of  a  material  provides  a  measure  of  how  fast  the  heat  transport  is  taking 
place  in  a  diffusive  manner. 

In  order  to  gain  further  understanding  regarding  behavior  and  dynamics  of  heat  diffusion,  we 
performed  thermal  relaxation  simulations  on  an  elongated  slab  with  the  length  of  -34 
nanometers.  Firstly,  the  entire  slab  was  annealed  and  equilibrated  at  very  low  temperature  to 
remove  any  history  of  thermal  fluctuations.  Thereafter,  the  20  A  region  in  the  centre  of  the  slab 
was  heated  to  desired  temperature  To.  Once  the  central  region  was  equilibrated  at  T0,  the 
thennostat  was  switched  off  and  constant  energy  simulations  were  performed  to  simulate  the  one 
dimensional  heat  flow  as  a  function  of  time.  The  simulations  were  performed  for  desired 
temperatures,  T0  of  100,  200,  300,  400,  500  and  600K.  Two  of  such  heat  flow  plots  are  shown  in 
Figure  26  and  27  across  the  whole  slab  as  a  function  of  time  for  first  500  picoseconds.  The  figure 
shows  that  the  temperature  of  the  central  region  decays  quite  sharply  while  the  temperature  of 
rest  of  the  slab  builds  up  over  time. 

The  temperature  decay  for  the  central  slab  is  better  depicted  in  Figure  28.  In  this  figure,  the 
decay  has  been  normalized  in  order  to  differentiate  the  relative  behavior  of  the  central  region.  It 
is  clear  from  the  figure  that  higher  initial  temperature  difference  leads  to  sharper  decay  of  central 
slab  temperature.  In  order  to  understand  the  decay  further,  the  plot  was  then  fitted  to  stretched 
exponential  curve  with  characteristic  relaxation  time,  x.  The  corresponding  stretching  factor  /? 
and  relaxation  time  r  are  shown  in  Table  6.  The  non-unity  values  of  y 3  suggest  that  there  exists  a 
distribution  of  relaxation  times  which  could  be  associated  with  the  low  frequency  vibrational 
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modes  necessary  for  heat  diffusion.  Although,  /?  increases  with  temperature,  its  non-unity  value 
suggests  that  heat  transport  is  indeed  a  phenomenon  associated  with  a  distribution  of  vibrational 
modes,  irrespective  of  temperature,  and  hence  a  distribution  of  relaxation  times.  The  table  also 
shows  the  x  decreases  as  the  temperature  difference  increases.  Such  a  decrease  could  suggest  the 
shift  of  relaxation  times  towards  lower  values,  associated  with  more  vibrant  motions  at  high 
temperatures.  It  is  also  interesting  to  see  that  the  distribution  behavior  of  relaxation  times  is 
much  more  sensitive  at  low  temperatures  (<  Tg)  than  high  temperature  (  >  Tg). 

Given  the  diffusivity  of  the  system,  it  is  possible  to  correlate  this  relaxation  time  x  with  a  length 
scale.  This  length  scale  should  provide  a  measure  of  how  far  the  thermal  transport  has  been 
reached  before  the  local  equilibrium  has  occurred,  as  associated  with  x.  In  order  to  observe  that 
length  scale  and  its  temperature  dependence  for  different  studied  T0,  we  used  one-dimensional 
analog  of  diffusivity  equation  shown  below. 

(/2)V2  =  V2Z)r 

To  observe  qualitative  trends,  thermal  diffusivity  was  assumed  to  be  independent  of  temperature. 
The  results  for  all  studied  cases  are  presented  in  Table  7.  The  table  shows  that  this  length  scale  is 
of  the  order  of  few  nanometers  and  decreases  with  increasing  temperature  difference.  Such  a 
decrease  could  be  attributed  to  more  number  of  vibrational  modes,  higher  probability  of 
vibrational  scattering  which  lead  to  low  values  of  relaxation  time  and  hence  low  /. 

The  temperature  dependence  of  this  length  scale  suggests  that  this  length  scale  should  be 
associated  with  mean  free  length  of  vibrational  modes  in  solids,  i.e.,  epoxy  networks  in  our  case. 
In  order  to  estimate  the  mean  free  length  from  experimental  values,  we  used  well  known  relation 
which  relates  diffusivity  to  the  mean  free  length  as  shown  below 


Here,  E,  p,  v  and  D  represent  elastic  modulus,  density,  sound  velocity  and  thennal  diffusivity  of 
the  material,  respectively.  For  epoxy  networks,  using  experimental  values  of  E,  p  and  D  as  ~2 
GPa,  1.17gm/cm  and  0.157  mm/s,  the  mean  free  path  /,  was  estimated  to  be  ~8  nanometers  at 
room  temperature.  It  is  interesting  to  observe  that  mean  free  length  and  observed  length  scale 
based  on  relaxation  time  and  diffusivity  calculations  are  at  same  order  of  magnitude.  It  should 
not  come  as  a  surprise  that  they  indeed  be  closely  associated.  The  former  length  scale  is 
associated  with  the  collision  of  low  frequency  vibrational  modes,  necessary  for  heat  diffusion 
while  the  later  provides  a  length  scale  of  heat  movement  in  order  to  drive  local  thermal 
equilibrium. 
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Figure  26.  Temperature  profile  of  the  entire  slab  as  a  function  of  time.  T0  =  600K. 
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Figure  28.  Normalized  temperature  decay  profile  of  the  central  region  for  different 
starting  T0.  Studied  temperatures  T0  were  100,  200,  300,  400,  500  and  600  K.  The  arrow 
points  towards  the  direction  of  increasing  temperature.  Inset  shows  the  stretched 
exponential  fitting  according  to  the  parameters  mentioned  in  Table  7. 


Table  6:  Stretched  exponential  fit6  parameters  and  calculated  length  scale 


Temperature 

To 

Stretching 

Coefficient 

P 

Characteristic 
Relaxation  Time 

T 

(/2f 

100 

0.39 

107.7 

57.5 

200 

0.41 

73.9 

47.7 

300 

0.43 

59.7 

42.8 

400 

0.45 

50.2 

39.3 

500 

0.44 

43.6 

36.6 

600 

0.45 

38.0 

34.1 

6  Fitted  to  the  equation  T(t)  =  e 
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Thermal  Conductivity  Simulations  for  CNT/Epoxy  Interface 


The  heat  transport  across  CNT/epoxy  network  interface  was  also  studied  using  molecular 
dynamics  simulations.  The  CNT/epoxy  network  system  which  was  build  to  perform  such 
simulations  is  shown  in  Figure  29  and  30.  For  these  simulations,  the  whole  box  was  divided  into 
12  shells  of  3  Angstroms  each  surrounding  the  CNT.  In  these  simulations,  the  CNT  at  the  centre 
of  the  slab  was  heated  continuously  while  same  amount  of  heat  was  taken  about  from  the 
outennost  shell.  After  a  significant  period  of  time  (-200  ps),  a  thermal  equilibrium  was  attained 
so  that  the  temperature  of  the  CNT  and  each  of  the  shells  was  not  changing  with  time.  The 
temperature  of  CNT  and  each  of  the  shells  was  then  calculated  from  velocities  and  was  block 
averaged  using  statistical  methods.  The  temperature  is  plotted  in  Figure  31.  It  is  very  clear  from 
the  figure  that  indeed  there  exists  a  huge  temperature  gradient  across  the  boundary  of  the 
CNT/epoxy  interface.  The  temperature  gradient  is  within  one  layer  of  epoxy  with  respect  to  CNT 
is  about  10  times  higher  than  temperature  gradient  across  the  whole  box  (  ~  40  Angs),  which 
suggest  the  dT/dx  at  the  interface  is  about  2  order  of  magnitude  higher  with  respect  to  epoxy 
materials.  Such  a  simulation  suggests  a  need  of  alternative  approaches  to  create  new  pathways 
for  the  heat  transport  within  epoxy  composite  materials. 


Figure  29:  Top  view  of  the  CNT/Epoxy  network  system. 
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Figure  30:  Slant  view  of  the  CNT/Epoxy  system 
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Figure  31:  Temperature  profile  across  CNT/epoxy  interface  and  within  epoxy  matrix. 
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